##### The following script consists of three parts:
### I. Mixed-effects models and robustness checks. This part replicates all models reported in the Supplementary material (Table A4).
### II. Figures. This code replicates all the figures reported in the main article and the Supplementary material.
### III. Appendix. This code replicates additional tables reported in the Supplementary material.



#####################################################
### I. Mixed-effects models and robustness checks ###
#####################################################

### MODELS

# The following script contains code for running the mixed-effect linear models using the lmer function from the lme4 package
# The models employ the full maximum likelihood estimator (REML = FALSE) with the Bound Optimization by Quadratic Approximation (bobyqa)
# bobyqa is used to achieve more stable convergence for complex models with several random effects estimated (cf Model 4)
# All models are weighted by the pspwght weight (post-stratified design weight)

library(lme4)
library(lmerTest)

# Model 0: empty model (random slopes and country-round and country levels)
M0_ld_pspwght_FML <- lmer(red_sup ~ (1 | cntry_round) + (1 | cntry), data = ess_4_ld_models, 
                          weights = pspwght, 
                          REML = FALSE,
                          control = lmerControl(optimizer = "bobyqa"))
summary(M0_ld_pspwght_FML)
AIC(M0_ld_pspwght_FML) # 948281.4
BIC(M0_ld_pspwght_FML) # 948324.1


### Model 1: only individual-level predictors
M1_ld_pspwght_FML <- lmer(red_sup ~ unemployed_CWC + employed_CWC + hinc_2_cop + hinc_3_dif + hinc_4_vdif + female + age + isced2 + isced3 + isced4 + isced56 + lrscale + factor(essround) + 
                            (1 | cntry_round) + (1 | cntry), data = ess_4_ld_models, weights = pspwght,
                          REML = FALSE,
                          control = lmerControl(optimizer = "bobyqa"))
summary(M1_ld_pspwght_FML)
AIC(M1_ld_pspwght_FML) # 924775.8
BIC(M1_ld_pspwght_FML) # 925032.3

# Model 2: adding three contextual-level predictors: unemployment rates, social expenditure and Gini coefficient of disposable income inequality
M2_ld_pspwght_FML <- lmer(red_sup ~ unemployed_CWC + employed_CWC + hinc_2_cop + hinc_3_dif + hinc_4_vdif + female + age + isced2 + isced3 + isced4 + isced56 + lrscale + factor(essround) + 
                            mean_HUR_q0 + dev_mean_HUR_q0 + mean_gini_disp_y0 + dev_mean_gini_disp_y0 + mean_SOCEXP_y0 + dev_mean_SOCEXP_y0 + 
                            (1 | cntry_round) + (1 | cntry), data = ess_4_ld_models, weights = pspwght,
                          REML = FALSE,
                          control = lmerControl(optimizer = "bobyqa"))
summary(M2_ld_pspwght_FML)
AIC(M2_ld_pspwght_FML) # 924762.9
BIC(M2_ld_pspwght_FML) # 925083.5
# Likelihood ratio test which compares the two models
anova(M1_ld_pspwght_FML, M2_ld_pspwght_FML)


# Model 3: adding random slope for WE of unemployment rates at the COUNTRY level
M3_ld_pspwght_FML <- lmer(red_sup ~ unemployed_CWC + employed_CWC + hinc_2_cop + hinc_3_dif + hinc_4_vdif + female + age + isced2 + isced3 + isced4 + isced56 + lrscale + factor(essround) + 
                            mean_HUR_q0 + dev_mean_HUR_q0 + mean_gini_disp_y0 + dev_mean_gini_disp_y0 + mean_SOCEXP_y0 + dev_mean_SOCEXP_y0 +
                            (1 | cntry_round) + (1 + dev_mean_HUR_q0 | cntry), data = ess_4_ld_models, weights = pspwght, 
                          REML = FALSE,
                          control = lmerControl(optimizer = "bobyqa"))
summary(M3_ld_pspwght_FML)
AIC(M3_ld_pspwght_FML) # 924743.1
BIC(M3_ld_pspwght_FML) # 925085.1
# Likelihood ratio test which compares the two models
anova(M2_ld_pspwght_FML, M3_ld_pspwght_FML)


# Model 4: adding the cross-level interaction between dev_mean_HUR_q0 and mean_SOCEXP_y0
# Random slope for WCE of unemployment rate are still present at the COUNTRY level
M4_ld_pspwght_FML <- lmer(red_sup ~ unemployed_CWC + employed_CWC + hinc_2_cop + hinc_3_dif + hinc_4_vdif + female + age + isced2 + isced3 + isced4 + isced56 + lrscale + factor(essround) +
                            mean_HUR_q0 + dev_mean_HUR_q0 + mean_gini_disp_y0 + dev_mean_gini_disp_y0 + mean_SOCEXP_y0 + dev_mean_SOCEXP_y0 +
                            dev_mean_HUR_q0:mean_SOCEXP_y0 + 
                            (1 | cntry_round) + (1 + dev_mean_HUR_q0 | cntry), data = ess_4_ld_models, weights = pspwght,
                          REML = FALSE,
                          control = lmerControl(optimizer = "bobyqa"))
summary(M4_ld_pspwght_FML)
AIC(M4_ld_pspwght_FML) # 924739.6
BIC(M4_ld_pspwght_FML) # 925092.3
# Likelihood ratio test which compares the two last models
anova(M3_ld_pspwght_FML, M4_ld_pspwght_FML)

library(sjPlot)
### Exporting the model results to html document
# The exported html document was used to create TABLE A4 in the Supplementary material (TABLE A4: Hierarchical linear models of redistribution support.)
tab_model(M0_ld_pspwght_FML, M1_ld_pspwght_FML, M2_ld_pspwght_FML, M3_ld_pspwght_FML, M4_ld_pspwght_FML, 
          digits = 3, digits.p = 4, digits.re = 4,
          show.se = TRUE, show.ci = FALSE, 
          show.aic = FALSE, 
          show.loglik = TRUE, show.dev = TRUE,
          string.est = "Estimate", string.se = "SE", 
          p.threshold = c(0.1, 0.05, 0.01),
          p.style = "stars")

# Comparing AICs of all five estimated models: the AICs gradually decrease as the model becomes more complex (the final model 4 providing the best fit to the data) 
rank(c(AIC(M0_ld_pspwght_FML), AIC(M1_ld_pspwght_FML), AIC(M2_ld_pspwght_FML), AIC(M3_ld_pspwght_FML), AIC(M4_ld_pspwght_FML)))



### ROBUSTNESS CHECKS
# These checks are run on the final Model 4: M4_ld_pspwght_FML

# (1) Excluding COUNTRIES (one at a time) from the model and estimating the final model j times on j-1 countries
# obtaining vector of countries on which the model is computed (from the model object: M4_ld_pspwght_FML)
countries_in_model <- levels(M4_ld_pspwght_FML@frame$cntry)
length(countries_in_model)

# creating a data frame where coefficients will be recorded: 3 coefficients and associated quantities are stored
robust_country <- data.frame(matrix(ncol = 15, nrow = length(countries_in_model)))
rownames(robust_country) <- countries_in_model
colnames(robust_country) <- c("mean_HUR_q0_Est", "mean_HUR_q0_SE", "mean_HUR_q0_df", "mean_HUR_q0_t", "mean_HUR_q0_p",
                              "dev_mean_HUR_q0_Est", "dev_mean_HUR_q0_SE", "dev_mean_HUR_q0_df", "dev_mean_HUR_q0_t", "dev_mean_HUR_q0_p",
                              "dev_mean_HUR_q0:mean_SOCEXP_y0_Est", "dev_mean_HUR_q0:mean_SOCEXP_y0_SE", "dev_mean_HUR_q0:mean_SOCEXP_y0_df", 
                              "dev_mean_HUR_q0:mean_SOCEXP_y0_y0_t", "dev_mean_HUR_q0:mean_SOCEXP_y0_p")

# assigning the coeffiecient estimates, their standard errors, degrees of freedom, t-values and p-vlues to the data frame
# !!! running the following script takes approximately 2 hours
for(i in 1:length(countries_in_model)){
  model <- lmer(red_sup ~ unemployed_CWC + employed_CWC + hinc_2_cop + hinc_3_dif + hinc_4_vdif + female + age + isced2 + isced3 + isced4 + isced56 + lrscale + factor(essround) + 
                  mean_HUR_q0 + dev_mean_HUR_q0 + mean_gini_disp_y0 + dev_mean_gini_disp_y0 + mean_SOCEXP_y0 + dev_mean_SOCEXP_y0 +
                  dev_mean_HUR_q0:mean_SOCEXP_y0 + 
                  (1 | cntry_round) + (1 + dev_mean_HUR_q0 | cntry), 
                data = ess_4_ld_models, subset = cntry != countries_in_model[i], weights = pspwght,
                REML = FALSE,
                control = lmerControl(optimizer = "bobyqa"))
  summary_model <- summary(model)
  coefficients <- summary_model$coefficients
  robust_country[i, 1:5] <- coefficients["mean_HUR_q0", ]
  robust_country[i, 6:10] <- coefficients["dev_mean_HUR_q0", ]
  robust_country[i, 11:15] <- coefficients["dev_mean_HUR_q0:mean_SOCEXP_y0", ]
  print(countries_in_model[i])
  print(summary(model))
}

# Saving the resulting quantities
write.csv(x = robust_country, "robust_country.csv")


##### (2) Excluding ESS-ROUNDS (one at a time) from the model and estimating the final model j times on j-1 ESS rounds
# creating vector of ESS rounds on which the final model is computed 
rounds_in_model <- seq(1:9)

# creating a data frame where coefficients and related quantities will be recorded
robust_round <- data.frame(matrix(ncol = 15, nrow = length(rounds_in_model)))
rownames(robust_round) <- rounds_in_model
colnames(robust_round) <- c("mean_HUR_q0_Est", "mean_HUR_q0_SE", "mean_HUR_q0_df", "mean_HUR_q0_t", "mean_HUR_q0_p",
                            "dev_mean_HUR_q0_Est", "dev_mean_HUR_q0_SE", "dev_mean_HUR_q0_df", "dev_mean_HUR_q0_t", "dev_mean_HUR_q0_p",
                            "dev_mean_HUR_q0:mean_SOCEXP_y0_Est", "dev_mean_HUR_q0:mean_SOCEXP_y0_SE", "dev_mean_HUR_q0:mean_SOCEXP_y0_df", 
                            "dev_mean_HUR_q0:mean_SOCEXP_y0_y0_t", "dev_mean_HUR_q0:mean_SOCEXP_y0_p")

# assigning the coeffiecient estimates and four related quantities to the data frame
# !!! running the following script takes approximately 40 minutes
for(i in 1:length(rounds_in_model)){
  model <- lmer(red_sup ~ unemployed_CWC + employed_CWC + hinc_2_cop + hinc_3_dif + hinc_4_vdif + female + age + isced2 + isced3 + isced4 + isced56 + lrscale + factor(essround) + 
                  mean_HUR_q0 + dev_mean_HUR_q0 + mean_gini_disp_y0 + dev_mean_gini_disp_y0 + mean_SOCEXP_y0 + dev_mean_SOCEXP_y0 +
                  dev_mean_HUR_q0:mean_SOCEXP_y0 + 
                  (1 | cntry_round) + (1 + dev_mean_HUR_q0 | cntry), 
                data = ess_4_ld_models, subset = essround != rounds_in_model[i], weights = pspwght,
                control = lmerControl(optimizer = "bobyqa"))
  summary_model <- summary(model)
  coefficients <- summary_model$coefficients
  robust_round[i, 1:5] <- coefficients["mean_HUR_q0", ]
  robust_round[i, 6:10] <- coefficients["dev_mean_HUR_q0", ]
  robust_round[i, 11:15] <- coefficients["dev_mean_HUR_q0:mean_SOCEXP_y0", ]
  print(rounds_in_model[i])
  print(summary(model))
}

# Saving the resulting quantities
write.csv(x = robust_round, "robust_essround.csv")





###################
### II. Figures ###
###################

##### Figure 1: Support for income redistribution and unemployment rates in Europe.

### Aggregating the data for the country-round time series figure (weighted average redistribution support is computed)
ess_cntry_round_figure <- ess_4_ld_models[, .(wght_mean_support = weighted.mean(x = red_sup, w = pspwght),
                                              HUR = mean(HUR_q0), dev_mean_HUR_q0 = mean(dev_mean_HUR_q0)), 
                                          by = c("cntry", "essround")]

countries <- unique(ess_cntry_round_figure$cntry)

### a. COUNTRY-ROUND SUPPORT FOR REDISTRIBUTION
# creating an empty matrix to construct the time series plots
full_matrix_red_sup <- matrix(data = rep(NA, 27*9), nrow = 27, ncol = 9)
rownames(full_matrix_red_sup) <- countries
colnames(full_matrix_red_sup) <- c("Round 1", "Round 2", "Round 3", "Round 4", "Round 5", 
                                   "Round 6", "Round 7", "Round 8", "Round 9")

# assignment of available aggregate country-round support for redistribution
for(i in 1:nrow(ess_cntry_round_figure)){
  country_i <- ess_cntry_round_figure$cntry[i]
  round_i <- as.numeric(ess_cntry_round_figure$essround[i])
  
  row_full <- as.numeric(which(rownames(full_matrix_red_sup) == country_i))
  column_full <- round_i
  
  full_matrix_red_sup[row_full, column_full] <- ess_cntry_round_figure$wght_mean_support[i]
}
# computing country-mean redistribution support
rowMeans(full_matrix_red_sup, na.rm = TRUE)

### b. COUNTRY-ROUND HARMONISED UNEMPLOYMENT RATE
# creating an empty matrix to construct the time series plots
full_matrix_HUR <- matrix(data = rep(NA, 27*9), nrow = 27, ncol = 9)
rownames(full_matrix_HUR) <- countries
colnames(full_matrix_HUR) <- c("Round 1", "Round 2", "Round 3", "Round 4", "Round 5", 
                               "Round 6", "Round 7", "Round 8", "Round 9")

# assignment of available aggregate country-round support for redistribution
for(i in 1:nrow(ess_cntry_round_figure)){
  country_i <- ess_cntry_round_figure$cntry[i]
  round_i <- as.numeric(ess_cntry_round_figure$essround[i])
  
  row_full <- as.numeric(which(rownames(full_matrix_HUR) == country_i))
  column_full <- round_i
  
  full_matrix_HUR[row_full, column_full] <- ess_cntry_round_figure$HUR[i]
}

# range of average support for income redistribution across 202 country rounds
range(ess_cntry_round_figure$wght_mean_support)

# range of HUR across 202 country rounds
range(ess_cntry_round_figure$HUR)

# country labels used in the figure
full_names <- c("Austria (AT)", "Belgium (BE)", "Bulgaria (BG)", "Switzerland (CH)", "Cyprus (CY)", 
                "Czech Republic (CZ)", "Germany (DE)", "Denmark (DK)", "Estonia (EE)", "Spain (ES)", 
                "Finland (FI)", "France (FR)", "United Kingdom (GB)", "Greece (GR)", "Hungary (HU)", 
                "Ireland (IE)", "Israel (IL)", "Iceland (IS)", "Italy (IT)", "Lithuania (LT)", 
                "Netherlands (NL)", "Norway (NO)", "Poland (PL)", "Portugal (PT)", "Sweden (SE)", 
                "Slovenia (SI)", "Slovakia (SK)")

### FIGURE 1: Creating the time series plot itself
### Exporting the plot: (PDF Size = Device Size: 12 x 8.13)
### Exporting the plot: (JPEG Size = 1400 x 800)
mat <- matrix(c(1:27, 28, 28, 28),
              nrow = 6, ncol = 5,
              byrow = TRUE)
layout(mat = mat)

par(las = 1, mar = c(2.5, 2.75, 1.75, 2.5))
for(i in 1:nrow(full_matrix_red_sup)){
  plot(seq(from = 1, to = 9, by = 1), full_matrix_red_sup[i, ], ylim = c(1.75, 3.85), main = full_names[i], 
       xaxt = "n", pch = 19, cex = 1.3, type = "o", xlab = "", ylab = "", lty = 3, lwd = 1.5, las = 1)
  axis(1, at = 1:9, labels = c(2002,2004, 2006, 2008, 2010, 2012, 2014, 2016, 2018))
  correlation <- round((cor(full_matrix_red_sup[i, ], full_matrix_HUR[i, ], use = "complete.obs")), digits = 2)
  text(x = 2.1, y = 3.7, labels = bquote(r == .(correlation)), cex = 1.1)
  par(new = TRUE)
  plot(seq(from = 1, to = 9, by = 1), full_matrix_HUR[i, ], ylim = c(0, 28), xaxt = "n", yaxt="n", pch = 15, 
       cex = 1.3, type = "o", xlab = "", ylab = "", lty = 1, col = "gray")
  axis(side = 4, las = 1, col.axis = "gray", col.ticks = "gray")
}

plot.new()
legend("bottomright", pch = c(19, 15), col = c("black", "gray"), lty = c(3, 1), 
       legend = c("Redistribution support (vertical axes on the left)", "Unemployment rate (vertical axes on the right)"), 
       bty = "n", cex = 1.5, pt.cex = c(1.5, 1.5), text.col = c("black", "gray"), y.intersp = 1.25)



##### Figure 2: Longitudinal and cross-sectional relationships between redistribution support and unemployment rates.

# Aggregated data for the COUNTRY-ROUND effect figure (WE)
par(mfrow = c(2, 1), mar = c(4.25, 4.5, 2.2, 0.5))
### Exporting the plot: (PDF Size = Device Size: 7.50 x 6.50)
### Exporting the plot: (JPEG Size = 720 x 586)

# Plot to the article
plot(ess_cntry_round_figure$dev_mean_HUR_q0, ess_cntry_round_figure$wght_mean_support, pch = 19, cex = 0.75,
     ylim = c(1.75, 4), xlim = c(-10, 10), 
     xlab = expression(paste("Deviation from the country−mean harmonised unemployment rate", ~ (HUR[tjCWC]))), 
     ylab = expression(paste("Redistribution support: country−round", ~~~ (bar(y)[tj]))), cex.lab = 0.85,
     las = 1, main = "Longitudinal relationship")
abline(lm(ess_cntry_round_figure$wght_mean_support ~ ess_cntry_round_figure$dev_mean_HUR_q0), lwd = 2.5, lty = 2, col = "grey")
cor(ess_cntry_round_figure$dev_mean_HUR_q0, ess_cntry_round_figure$wght_mean_support, use = "complete.obs")
text(x = 8, y = 1.9, labels = expression(paste("r = 0.06, ", ~ n[tj] == 202)), cex = 1.25)
mtext(expression(bar(y)[tj] == 2.857 + 0.007*HUR[tjCWC]), side = 3, line = -2,  col = "gray", cex = 1.1)
lm(ess_cntry_round_figure$wght_mean_support ~ ess_cntry_round_figure$dev_mean_HUR_q0)

# Aggregated data for the COUNTRY effects figures (BE)
# computing the weighted mean-country support for redistribution ("mean_support_w") - each weighted mean is based on 4-9 country-round means (i.e. this procedure aggregates the country-round means)
ess_cntry_figure <- ess_cntry_round_figure[, .(mean_support_w = mean(wght_mean_support), mean_HUR = mean(HUR)), by = .(cntry)]

# Plot to the article
plot(ess_cntry_figure$mean_HUR, ess_cntry_figure$mean_support_w, pch = 19, cex = 1.1,
     ylim = c(1.8, 4), xlim = c(0, 20), 
     xlab = expression(paste("Average harmonised unemployment rate", ~~ (bar(HUR)[j]))), 
     ylab = expression(paste("Redistribution support: country", ~~ (bar(y)[j]))), las = 1, cex.lab = 0.85,
     main = "Cross−sectional relationship")
abline(lm(ess_cntry_figure$mean_support_w ~ ess_cntry_figure$mean_HUR), lwd = 2, lty = 2, col = "gray")
text(ess_cntry_figure$mean_HUR, ess_cntry_figure$mean_support_w, labels = ess_cntry_figure$cntry, pos = 3, cex = 0.7)
cor(ess_cntry_figure$mean_HUR, ess_cntry_figure$mean_support_w)
text(18.25, 2, labels = expression(paste("r = 0.53, ", ~ n[j] == 27)), cex = 1.25)
mtext(expression(bar(y)[j] == 2.429 +0.058*bar(HUR)[j]), side = 3, line = -2.25,  col = "gray", cex = 1.1)
lm(ess_cntry_figure$mean_support_w ~ ess_cntry_figure$mean_HUR)



##### Figure 3: FIGURE 3. Point estimates and confidence intervals for the effects of contextual variables.
# Storing the respective coefficient estimates and their standard errors from models 2, 3 and 4 into vectors
coefs_2 <- fixef(M2_ld_pspwght_FML)[c("mean_HUR_q0", "dev_mean_HUR_q0", "mean_gini_disp_y0", "dev_mean_gini_disp_y0", "mean_SOCEXP_y0", "dev_mean_SOCEXP_y0")]
se_2 <- sqrt(diag(vcov(M2_ld_pspwght_FML)))[c("mean_HUR_q0", "dev_mean_HUR_q0", "mean_gini_disp_y0", "dev_mean_gini_disp_y0", "mean_SOCEXP_y0", "dev_mean_SOCEXP_y0")]
coefs_3 <- fixef(M3_ld_pspwght_FML)[c("mean_HUR_q0", "dev_mean_HUR_q0", "mean_gini_disp_y0", "dev_mean_gini_disp_y0", "mean_SOCEXP_y0", "dev_mean_SOCEXP_y0")]
se_3 <- sqrt(diag(vcov(M3_ld_pspwght_FML)))[c("mean_HUR_q0", "dev_mean_HUR_q0", "mean_gini_disp_y0", "dev_mean_gini_disp_y0", "mean_SOCEXP_y0", "dev_mean_SOCEXP_y0")]
coefs_4 <- fixef(M4_ld_pspwght_FML)[c("mean_HUR_q0", "dev_mean_HUR_q0", "mean_gini_disp_y0", "dev_mean_gini_disp_y0", "mean_SOCEXP_y0", "dev_mean_SOCEXP_y0", "dev_mean_HUR_q0:mean_SOCEXP_y0")]
se_4 <- sqrt(diag(vcov(M4_ld_pspwght_FML)))[c("mean_HUR_q0", "dev_mean_HUR_q0", "mean_gini_disp_y0", "dev_mean_gini_disp_y0", "mean_SOCEXP_y0", "dev_mean_SOCEXP_y0", "dev_mean_HUR_q0:mean_SOCEXP_y0")]

# Storing variable labels into a vector
labely <- c("Unemployment rate (BE)", 
            "Unemployment rate (WE)", 
            "Gini coefficient of disposable income (BE)", 
            "Gini coefficient of disposable income (WE)", 
            "Public social expenditure (BE)", 
            "Public social expenditure (WE)",
            "Unemp. rate (WE) * Public social exp. (BE)")

### Exporting the plot: (PDF Size = Device Size: 8.5 x 6.5)
### Exporting the plot: (JPEG Size = 820 x 590)

par(mfrow = c(1, 1), mar = c(4.5, 16.75, 3, 0.25))

# location of points on the vertical axis            
y <- seq(from = 7*1.5, to = 2, by = -1.5)
y_int <- seq(from = 7*1.5, to = 0.5, by = -1.5)

## Model 2: light gray colour
plot(x = coefs_2, y = y+0.3, ylim = c(1, 11), pch = 18, cex = 1.4, xlim = c(-0.05, 0.12), 
     ylab = "", yaxt = "n", xlab = "Contextual effects: coefficient estimates",
     col = "gray90")
axis(side = 3)
abline(v = 0, lty = 2)
text(par("usr")[1], y_int, labels = labely, pos = 2, xpd = TRUE, cex = 1)
for(i in 1:6){
  lines(x = c(coefs_2[i] - qnorm(0.975)*se_2[i], coefs_2[i] + qnorm(0.975)*se_2[i]),
        y = c(y[i]+0.3, y[i]+0.3), lwd = 1,
        col = "gray90")
}
for(i in 1:6){
  lines(x = c(coefs_2[i] - qnorm(0.95)*se_2[i], coefs_2[i] + qnorm(0.95)*se_2[i]),
        y = c(y[i]+0.3, y[i]+0.3), lwd = 2.5,
        col = "gray90")
}

## Model 3: dark gray colour
points(x = coefs_3, y = y, ylim = c(0, 7), pch = 15, cex = 1.1, col = "gray40")
for(i in 1:6){
  lines(x = c(coefs_3[i] - qnorm(0.975)*se_3[i], coefs_3[i] + qnorm(0.975)*se_3[i]),
        y = c(y[i], y[i]), lwd = 1,
        col = "gray40")
}
for(i in 1:6){
  lines(x = c(coefs_3[i] - qnorm(0.95)*se_3[i], coefs_3[i] + qnorm(0.95)*se_3[i]),
        y = c(y[i], y[i]), lwd = 2.5,
        col = "gray40")
}


## Model 4: black colour
points(x = coefs_4, y = y_int-0.3, ylim = c(0, 7), pch = 16, cex = 1.2)
for(i in 1:6){
  lines(x = c(coefs_4[i] - qnorm(0.975)*se_4[i], coefs_4[i] + qnorm(0.975)*se_4[i]),
        y = c(y_int[i]-0.3, y_int[i]-0.3), lwd = 1)
}
for(i in 1:6){
  lines(x = c(coefs_4[i] - qnorm(0.95)*se_4[i], coefs_4[i] + qnorm(0.95)*se_4[i]),
        y = c(y_int[i]-0.3, y_int[i]-0.3), lwd = 2.5)
}

# Adding legend to the figure
legend("bottomright",
       cex = 1,
       bty = "o",
       legend = c("Model 2", "Model 3", "Model 4"),
       col = c("gray90", "gray40", "black"),
       pch = c(18, 15, 16),
       pt.cex = c(1.4, 1.1, 1.2),
       lty = c(1, 1, 1),
       lwd = c(2.5, 2.5, 2.5))



##### Figure 4: Marginal effect of cyclical unemployment and public social expenditure

### Plotting the cross-level interaction (i.e. the marginal effect of cyclical unemployment on aggregate redistribution support at different welfare state sizes)
library(interplot)
### Exporting the plot: (PDF Size = Device Size: 8.5 x 7)
### Exporting the plot: (JPEG Size = 820 x 590)
interplot(m = M4_ld_pspwght_FML, 
          var1 = "dev_mean_HUR_q0", var2 = "mean_SOCEXP_y0", 
          ci = 0.95, hist = FALSE, adjCI = TRUE, sims = 10000) + 
  geom_line(linewidth = 1) + 
  xlab("Country-mean public social expenditure (percentage of GDP)") +
  ylab("Marginal effect: longitudinal effect (WE) of the unemployment rate") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black")) +
  geom_hline(yintercept = 0, linetype = "dotted", linewidth = 1) +
  ### The marginal effect of WE predicted by the model M4_ld_pspwght_FML for each country:
  # coefficient point estimate + cross-level interaction * public social expenditure + random effect estimate for each country
  # the following script implements this logic when computing the y value in the aes argument
  geom_point(data = ess_4_context_means, 
             aes(x = mean_SOCEXP_y0, y = fixef(M4_ld_pspwght_FML)["dev_mean_HUR_q0"] + fixef(M4_ld_pspwght_FML)["dev_mean_HUR_q0:mean_SOCEXP_y0"] * mean_SOCEXP_y0 + ranef(M4_ld_pspwght_FML)$cntry[, 2])) +
  geom_text(data = ess_4_context_means, aes(x = mean_SOCEXP_y0, y = fixef(M4_ld_pspwght_FML)["dev_mean_HUR_q0"] + fixef(M4_ld_pspwght_FML)["dev_mean_HUR_q0:mean_SOCEXP_y0"] * mean_SOCEXP_y0 + ranef(M4_ld_pspwght_FML)$cntry[, 2], label = cntry), 
            hjust = -0.35, vjust = 0.2, size = 3.25) + 
  ylim(-0.0723, 0.04725) +
  xlim(15.3, 30.7)



##### Figure A1: Jackknife estimates of two key coefficients (countries excluded).

### plotting the coefficients of Between Effect of unemployment rate
# Device Size: 10.00 x 7.00 (Portrait)
range(robust_country[, "mean_HUR_q0_Est"])
par(las = 1, mfrow = c(2, 1), mar = c(4.25, 4.25, 1, 0.5))
plot(robust_country[, "mean_HUR_q0_Est"], ylim = c(-0.02, 0.11), pch = 19, 
     ylab = "Unemployment rate (BE)", xlab = "Country excluded from the model", 
     xaxt = "n", cex = 1.1, main = NA)
axis(1, at = 1:27, labels = countries_in_model, cex.axis = 0.75)
abline(a = fixef(M4_ld_pspwght_FML)[22], b = 0, lty = 2, lwd = 2)
abline(a = 0, b = 0, lty = 1, lwd = 2)
legend("topleft", legend = "Coefficient estimate from the final model (i.e. Model 4 estimated on all 27 countries)", lty = 2, bty = "n", lwd = 2)
rect(xleft = 0, ybottom = -0.025, xright = 29, ytop = 0.0, density = NULL, angle = 45,
     col= rgb(0.6, 0.6, 0.6, alpha = 0.2), border = NA, lty = par("lty"), lwd = par("lwd"))

# computing the lower bound ($lb) and upper bound ($ub) of the 95% confidence intervals
robust_country$ub <- robust_country$mean_HUR_q0_Est + qt(p = 0.975, df = robust_country$mean_HUR_q0_df)*robust_country$mean_HUR_q0_SE
robust_country$lb <- robust_country$mean_HUR_q0_Est - qt(p = 0.975, df = robust_country$mean_HUR_q0_df)*robust_country$mean_HUR_q0_SE
# adding the 95% confidence intervals of coefficient estimates to the graph
for(i in 1:length(countries_in_model)){
  lines(c(i, i), c(robust_country[i, "lb"], robust_country[i, "ub"]))
}



##### Figure A2: Jackknife estimates of two key coefficients (ESS rounds excluded).

### plotting the coefficients of cross-level interaction
# Device Size: 10.00 x 7.00 (Portrait)
range(robust_country[, "dev_mean_HUR_q0:mean_SOCEXP_y0_Est"])
par(las = 1, mar = c(4.25, 4.25, 1, 0.5))
plot(robust_country[, "dev_mean_HUR_q0:mean_SOCEXP_y0_Est"], ylim = c(-0.01, 0.005), pch = 19, 
     ylab = "Unemployment rate (WE) * Social expenditure (BE)", 
     xlab = "Country excluded from the model", 
     xaxt = "n", cex.lab = 0.725, cex.axis = 0.75)
axis(1, at = 1:27, labels = countries_in_model, cex.axis = 0.75)
abline(a = fixef(M4_ld_pspwght_FML)[28], b = 0, lty = 2, lwd = 2)
abline(a = 0, b = 0, lty = 1, lwd = 2)
legend("bottomleft", legend = "Estimate of the cross-level interaction from the final model (i.e. Model 4 estimated on all 27 countries)", lty = 2, bty = "n", lwd = 2)
rect(xleft = 0, ybottom = 0, xright = 29, ytop = 0.25, density = NULL, angle = 45,
     col= rgb(0.6, 0.6, 0.6, alpha = 0.2), border = NA, lty = par("lty"), lwd = par("lwd"))

# computing the lower bound ($lb) and upper bound ($ub) of the 95% confidence intervals
robust_country$ub <- robust_country$`dev_mean_HUR_q0:mean_SOCEXP_y0_Est` - qt(p = 0.975, df = robust_country$`dev_mean_HUR_q0:mean_SOCEXP_y0_df`)*robust_country$`dev_mean_HUR_q0:mean_SOCEXP_y0_SE`
robust_country$lb <- robust_country$`dev_mean_HUR_q0:mean_SOCEXP_y0_Est` + qt(p = 0.975, df = robust_country$`dev_mean_HUR_q0:mean_SOCEXP_y0_df`)*robust_country$`dev_mean_HUR_q0:mean_SOCEXP_y0_SE`
# adding the 95% confidence intervals of coefficient estimates to the graph
for(i in 1:length(countries_in_model)){
  lines(c(i, i), c(robust_country[i, "lb"], robust_country[i, "ub"]))
}





#####################
### III. Appendix ###
#####################

##### Table A1: Sample sizes in analysed country rounds
table_A1 <- table(ess_4_ld_models$cntry, ess_4_ld_models$essround)
dim(ess_4_ld_models)
write.csv(x = table_A1, file = "Table_A1_appendix.csv")


##### Table A2: Country means of individual-level explanatory variables
describe_individual <- ess_4_ld_models[, .(
  unemployed = round(mean(unemployed), digits = 3), 
  employed = round(mean(employed), digits = 3),
  hinc_2_cop = round(mean(hinc_2_cop), digits = 3), 
  hinc_3_dif = round(mean(hinc_3_dif), digits = 3), 
  hinc_4_vdif = round(mean(hinc_4_vdif), digits = 3),
  female = round(mean(female), digits = 3), 
  age = round(mean(age), digits = 2), 
  isced2 = round(mean(isced2), digits = 3), 
  isced3 = round(mean(isced3), digits = 3), 
  isced4 = round(mean(isced4), digits = 3), 
  isced56 = round(mean(isced56), digits = 3), 
  lrscale = round(mean(lrscale), digits = 2)
), by = "cntry"]

write.csv(x = describe_individual, file = "Table_A2_appendix.csv")


##### TABLE A3: Contextual data at the country-round (i.e. values of contextual variables at country-round level)
table_A3 <- ess_4_ld_models[, .(quarter = key_quartal[1], HUR = mean(HUR_q0), SOC_EXP = mean(SOCEXP_y0), GINI = mean(gini_disp_y0)), by = c("cntry", "essround")]
dim(table_A3)

write.csv(table_A3, file = "Table_A3_appendix.csv")


##### TABLE A5: Hierarchical linear models of redistribution support: including GDP per capita (in USD) as an explanatory variable.
### Estimating the models with two variables measuring GDP: mean_GDP_y0 and dev_mean_GDP_y0
# NOTE: Models 0 and 1 remain unchanged (since they do not contain any variables measured at higher levels). Hence, models 0 and 1 correspond to the Tabla A4.

# Model 2: adding four contextual-level predictors: unemployment rates, social expenditure, Gini coefficient of disposable income inequality and GDP per capita
M2_ld_pspwght_FML_gdp <- lmer(red_sup ~ unemployed_CWC + employed_CWC + hinc_2_cop + hinc_3_dif + hinc_4_vdif + female + age + isced2 + isced3 + isced4 + isced56 + lrscale + factor(essround) + 
                                mean_HUR_q0 + dev_mean_HUR_q0 + mean_gini_disp_y0 + dev_mean_gini_disp_y0 + mean_SOCEXP_y0 + dev_mean_SOCEXP_y0 + mean_GDP_y0 + dev_mean_GDP_y0 +
                                (1 | cntry_round) + (1 | cntry), data = ess_4_ld_models, weights = pspwght,
                              REML = FALSE,
                              control = lmerControl(optimizer = "bobyqa"))
summary(M2_ld_pspwght_FML_gdp)
AIC(M2_ld_pspwght_FML_gdp) # 924765
BIC(M2_ld_pspwght_FML_gdp) # 925107

# Model 3: adding random slope for WE of unemployment rates at the COUNTRY level
M3_ld_pspwght_FML_gdp <- lmer(red_sup ~ unemployed_CWC + employed_CWC + hinc_2_cop + hinc_3_dif + hinc_4_vdif + female + age + isced2 + isced3 + isced4 + isced56 + lrscale + factor(essround) + 
                                mean_HUR_q0 + dev_mean_HUR_q0 + mean_gini_disp_y0 + dev_mean_gini_disp_y0 + mean_SOCEXP_y0 + dev_mean_SOCEXP_y0 + mean_GDP_y0 + dev_mean_GDP_y0 +
                                (1 | cntry_round) + (1 + dev_mean_HUR_q0 | cntry), data = ess_4_ld_models, weights = pspwght, 
                              REML = FALSE,
                              control = lmerControl(optimizer = "bobyqa"))
summary(M3_ld_pspwght_FML_gdp)
AIC(M3_ld_pspwght_FML_gdp) # 924743.8
BIC(M3_ld_pspwght_FML_gdp) # 925107.1

# Model 4: adding the cross-level interaction between dev_mean_HUR_q0 and mean_SOCEXP_y0
# Random slope for WCE of unemployment rate are still present at the COUNTRY level
M4_ld_pspwght_FML_gdp <- lmer(red_sup ~ unemployed_CWC + employed_CWC + hinc_2_cop + hinc_3_dif + hinc_4_vdif + female + age + isced2 + isced3 + isced4 + isced56 + lrscale + factor(essround) +
                                mean_HUR_q0 + dev_mean_HUR_q0 + mean_gini_disp_y0 + dev_mean_gini_disp_y0 + mean_SOCEXP_y0 + dev_mean_SOCEXP_y0 + mean_GDP_y0 + dev_mean_GDP_y0 +
                                dev_mean_HUR_q0:mean_SOCEXP_y0 + 
                                (1 | cntry_round) + (1 + dev_mean_HUR_q0 | cntry), data = ess_4_ld_models, weights = pspwght,
                              REML = FALSE,
                              control = lmerControl(optimizer = "bobyqa"))
summary(M4_ld_pspwght_FML_gdp)
AIC(M4_ld_pspwght_FML_gdp) # 924740.4
BIC(M4_ld_pspwght_FML_gdp) # 925114.4

### Exporting the model results to html document
# The exported html document was used to create TABLE A5 in the Supplementary material (TABLE A5. Hierarchical linear models of redistribution support: including GDP per capita (in USD).)
tab_model(M0_ld_pspwght_FML, M1_ld_pspwght_FML, M2_ld_pspwght_FML_gdp, M3_ld_pspwght_FML_gdp, M4_ld_pspwght_FML_gdp, 
          digits = 3, digits.p = 4, digits.re = 4,
          show.se = TRUE, show.ci = FALSE, 
          show.aic = FALSE, 
          show.loglik = TRUE, show.dev = TRUE,
          string.est = "Estimate", string.se = "SE", 
          p.threshold = c(0.1, 0.05, 0.01),
          p.style = "stars")
